#Plot of core and edge proportions (plus sds) with no parasite cost at k=32000 and cell.size =5

setwd("/Users/ben/Documents/Papers/Submitted/Rhabdias/Parasite lag model/Pathol Freq Dep Dyn Load")


y.error.bars<-function(x, y, ybar, lty="solid", ...){
	points(x, y, ...)
	arrows(x, y-ybar, x, y+ybar, code=3, angle=90, length=0.1, lty=lty, ...)
}

data<-read.table(file="Pathol Freq Depdt SI r15 Final.txt", header=T)
names(data)[3]<-"beta"
attach(data)

par(mfrow=c(2, 1)) #resize to 530 x 775
## First plot
core<-tapply(core.prop, list(surv.cost, beta), mean, na.rm=T)
edge<-tapply(edge.prop, list(surv.cost, beta), mean)
core.sd<-tapply(core.prop, list(surv.cost, beta), sd)
edge.sd<-tapply(edge.prop, list(surv.cost, beta), sd)
x<-as.numeric(colnames(core))

plot(c(min(x), max(x)), 0:1, ,type="n", bty="l", xlab="", ylab="Observed prevalence", main="Density dependent transmisison")
cols<-c("black", "grey40", "grey70")
for (i in 1:3){
	if (i==1) y.error.bars(x, core[i,], core.sd[i,], lty="solid", pch=19,cex=1.5, col=cols[i])
	else points(x, core[i,], pch=19, col=cols[i], cex=1.5)
	lines(x, core[i,], col=cols[i])
}
for (i in 1:3){
	if (i==1) y.error.bars(x, edge[i,], edge.sd[i,], lty="dashed", pch=21,cex=1.5, col=cols[i])
	else points(x, edge[i,], pch=21, col=cols[i], cex=1.5)
	lines(x, edge[i,], col=cols[i])
}

legend(0.12, 0.67, legend=c("Range core", "Range edge"), pch=c(19, 21),bty="n")
legend(0.12, 0.9, legend=c("=0", "=0.1", "=0.4"), lty="solid",col=cols, bty="n", title="Cost of parasitism")


## Second plot

summary<-tapply(x.diff, list(surv.cost, beta), mean)
summary[which(summary==Inf)]<-NA
summary.sd<-tapply(x.diff, list(surv.cost, beta), sd)
summary.sd[which(is.nan(summary.sd)==T)]<-NA

plot(c(min(x), max(x)), c(-4, 200), ,type="n", bty="l", yaxp=c(0, 200, 5), xlab="Transmission coefficient", ylab="Lag distance (arbitrary units)")
for (i in 1:3){
	if (i==1) y.error.bars(x, summary[i,], summary.sd[i,], lty="solid", pch=21,cex=1.5, col=cols[i])
	else points(x, summary[i,], pch=21, col=cols[i], cex=1.5)
	lines(x, summary[i,], col=cols[i])
}
legend(0.12, 200, legend=c("=0", "=0.1", "=0.4"), lty="solid",col=cols, bty="n", title="Cost of parasitism")
## Clean up

detach(data)
rm(data, core, core.sd, edge, edge.sd, x, y.error.bars, summary, summary.sd)

